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Numerical  simulation  is  a  fundamental  instrument  for  the  elaboration  and  assessment  of  a  strategic  utilization 
of  geothermal  energy.  It  can  be  used  for  the  evaluation  of  both  the  natural  (unperturbed)  state  and  the 
production  scenarios.  The  motivation  and  important  role  of  the  numerical  models  are  described  here  and 
deeply  illustrated  in  the  context  of  the  geothermal  energy  exploitation.  The  mathematical-physical  background 
is  also  briefly  illustrated,  together  with  all  the  practical  problems  of  modeling  and  implementation.  Particular 
attention  must  be  paid  to  the  boundary  conditions  and  thermophysical  parameters  assignment  and  calibration. 
The  reliability  of  the  model  must  be  accurately  evaluated,  in  order  to  prevent  common  failures  in  design  and 
running  of  the  energy  conversion  units  and  wells.  Several  case  studies  are  reviewed  and  discussed,  and  a  final 
discussion  is  presented.  The  limits  of  the  reservoir  modeling  and  simulation  are  also  outlined  in  a  general 
methodological  perspective  of  integrated  analysis.  The  scenarios  modeled  and  assessed  can  be  then  used  as 
practical  tools  for  the  sizing  and  optimization  of  the  power  unit  or  direct  heat  utilization. 
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1.  Introduction 

Geothermal  energy  is  considered  a  strategic  resource  in  many 
countries,  even  if  its  use  appears  to  be  often  marginal  in  the 
national  energy  systems.  Its  continuous  operating  mode  distin¬ 
guishes  geothermal  energy  from  the  other  renewable  sources, 
intermittent  or  stochastic.  Majority  of  the  geothermal  sources 
worldwide  are  medium-low  enthalpy  type  (water  dominant,  at 
temperature  lower  than  150  DC  and  pressure  below  15  bars). 
Stefansson,  [1]  estimated  that  more  than  70%  of  the  geothermal 
resources  available  in  the  world  are  water  dominated  fields,  at 
temperatures  under  150  °C.  The  binary  cycle  technology  with 
Organic  Rankine  Cycle  (ORC)  appears  to  be  the  most  efficient 
and  convenient  solution  for  such  a  kind  of  resource  [2],  Binary 
power  plants  are  now  objects  of  wide  attention  by  energy  markets, 
although  their  diffusion  is  still  made  difficult  by  a  lacking  technol¬ 
ogy  standardization  and  due  to  the  quite  high  specific  costs  [3,4], 
The  great  variability  of  the  resource  characteristics  worldwide  is 
one  of  the  possible  reasons.  The  proper  matching  between  the 
reservoir  capability  and  the  plants  parameters  (power  size,  extrac¬ 
tion/reinjection  rate)  is  a  critical  key  point.  In  power  plants  using 
dry-steam  (high  enthalpy)  geothermal  resources,  pressure  and 
temperature  reduction  can  be  compensated  by  an  increase  of  the 
mass  flow  rate.  In  case  of  binary  plants,  a  variation  of  the  resource 
properties  ( Tgeo ,  pgeo )  could  also  lead  to  a  fast  end  of  life  of  the 
plant.  The  first  and  most  important  activity  to  design  a  geothermal 
energy  plant  is  an  accurate  investigation  of  the  geothermal 
potential  assessment,  as  well  as  the  prediction  of  reservoir 
response  at  given  industrial  exploitation  configurations.  For  these 
reasons,  a  multidisciplinary  approach  to  the  problem  of  exploita¬ 
tion  of  geothermal  fields  (in  particular  at  medium-low  tempera¬ 
ture)  is  necessary.  Thermal  engineering,  geochemistry,  geophysics, 
and  reservoir  engineering  are  the  fields  involved  in  this  technique 
(Fig.  1).  The  authors  have  diffusely  discussed  this  topic  in  a  recent 
paper  [5], 

Numerical  simulation  is  a  fundamental  and  strongly  interacting 
instrument  for  plant  design  [6],  Different  approaches  are  here 
considered  with  reference  to  several  case  studies  of  geothermal  fields, 
which  are  reviewed  and  discussed.  The  perspectives  of  numerical 
simulation  of  geothermal  reservoirs  as  support  to  the  design  and 
sizing  of  geothermal  plants  are  also  outlined.  Simulation  can  be  very 
important  in  order  to  define  and  progressively  modify  the  manage¬ 
ment  strategy  of  the  geothermal  field.  Construction  of  the  numerical 
model  must  be  supported  by  a  detailed  knowledge  of  the  spatial 
distribution  of  the  properties  of  the  reservoir:  the  accuracy  in  the 
definition  of  the  dataset  is  fundamental  for  the  construction  of  the 
model.  The  model  is  then  enriched  by  including  the  database  of 
historical  data  collected  during  exploration. 

The  results  obtained  depend  a  lot  on  the  accuracy  level  of  the 
input  data.  The  model  will  be  much  more  accurate  if  as  much 
details  as  possible  are  known  about  the  geological  properties  of 
the  rocks  (effective  porosity,  density,  specific  heat,  permeability, 
thermal  conductivity),  thermophysical  properties  of  the  fluid 
(specific  heat,  density,  thermal  conductivity),  fractures  pattern 


/  \ 

Fig-1.  The  multidisciplinary  approach  proposed,  with  the  connections  between  the 
three  areas  involved. 


and  layout,  natural  recharge  of  fluid,  geothermal  boundary 
conditions. 

The  numerical  model  of  a  geothermal  reservoir  is  very  impor¬ 
tant  both  for  the  definition  of  the  geothermal  potential  assessment 
and  of  the  reinjection  strategy.  The  geothermal  potential  of  a 
particular  area  means  the  definition  of  temperature  ( Tgeo )  and 
pressure  (pgeo)  of  the  geothermal  fluid  and  also  of  the  maximum 
mass  flow  rate  (Mgeo)  that  can  be  extracted  maintaining  the 
thermal  properties  of  the  reservoir  and  of  the  geofluid  constant 
for  a  long  time.  Concerning  the  reinjection  strategy,  it  is  necessary 
to  take  into  account  the  circulation  model  of  the  fluid  in  the 
regional  area  considered  [7],  A  general  methodology  for  the 
reinjection  technologies  is  not  properly  available  in  literature, 
the  optimal  strategy  is  in  fact  site-dependant,  as  the  potential 
assessment  itself.  Interesting  discussions  on  this  particular  topic 
are  reported  by  SigurSsson  et  al.  [8],  and  recently  by  Kaya  et  al.  [9], 

The  main  task  of  potential  assessment  and  sustainable  plants 
design  is  the  optimization  and  enhancement  of  the  resource 
durability  (Vaccaro  [6]).  Interesting  contribution  on  the  definition 
and  evaluation  of  the  sustainability  and  renewability  of  the 
geothermal  energy  uses  are  available  in  the  recent  paper  by 
Hahnlein  et  al.  [10],  together  with  the  paper  by  Axelsson  [11], 
Particularly  in  case  of  innovative  geothermal  utilizations  (like  for 
example  EGS)  the  long-term  consequences  on  the  environment 
are  not  completely  known  yet.  The  same  argument  can  be  then 
referred  to  the  renewability  of  the  resource  itself,  which  is  directly 
dependent  on  the  type  and  rate  of  utilization.  The  renewability 
(and  sustainability)  reference  level  can  vary,  as  one  can  adjust  the 
energy  system  size  and  extraction  rate  according  to  an  acceptable 
durability  level  [11].  Also  in  case  of  direct  heat  utilization  (e.g. 
district  heating)  some  strategy  remarks  should  be  pointed  out.  In 
the  recent  paper  by  Fox  et  al.  [12]  the  renewable  capacity  of  deep 
systems  is  assessed  and  discussed  in  order  to  elaborate  a  rotating 
utilization  strategy. 

The  numerical  simulation  of  geothermal  reservoirs  is  a  well 
known  topic  and  has  already  been  an  object  of  investigations  and 
reviews  (e.g.  O'Sullivan  [13]).  Unfortunately,  till  now  the  use  of 
numerical  simulation  has  not  faced  any  direct  connection  with  the 
energy  systems  analysis.  A  proper  prediction  should  deal  with  the 
changes  of  the  different  parameters  in  response  to  given  mass  flow 
rate  extraction  and  reinjection  (corresponding  to  the  specific 
energy  strategy).  It  is  evident  that  a  key  role  is  assigned  here  to 
the  numerical  simulation  of  the  reservoir,  as  compared  to  other 
reservoir  engineering  aspects  (wells  siting,  fluid  losses,  tracer  test). 
This  ambitious  task  seems  to  be  not  of  specific  interest  in  most  of 
the  analyses  carried  out  in  the  past,  so  the  authors  would  like  to 
review  the  recent  developments  in  the  field  of  numerical  simula¬ 
tion  of  geothermal  fields,  focusing  their  attention  on  the  particular 
use  of  such  an  important  instrument  for  the  sustainable  design  of 
geothermal  plants. 


2.  Numerical  simulation  of  geothermal  reservoirs:  strategic 
role  for  the  design  of  geothermal  plants 

The  numerical  simulation  of  a  geothermal  reservoir  is  a  well 
known  field  of  research  in  the  literature  and  it  has  already  been  an 
object  of  accurate  review  analysis  and  methodological  overview 
[13-17], 

Two  main  goals  can  be  identified:  history  matching  and 
forecast  of  future  scenarios  (consequent  to  the  exploitation  of 
the  reservoir).  History  matching  is  usually  done  to  check  the 
reliability  of  a  model  and  evaluate  the  sustainability  level  in 
retrospect.  It  is  the  analysis  of  an  exploitation  history  according 
to  the  data  log  until  present  time  or  during  a  particular  time 
interval.  This  also  allows  to  check  the  numerical  model  in  a 
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“feedback  loop"  and  calibrate  or  adapt  it  to  what  is  reported  in  the 
history  data  (in  order  to  improve  future  scenarios  forecast).  This  is 
practically  done  by  assigning  temperature  and  mass  flow  rate  of 
both  the  geofluid  extracted  and  reinjected,  and  any  other  useful 
historical  data  recorded  that  can  be  translated  into  thermophysical 
parameters  or  boundary/initial  conditions.  Some  phenomena 
which  affect  the  behavior  of  field  utilization  could  not  be  put  into 
a  simulation  in  a  proper  way,  without  knowing  about  them  from 
the  history  (natural  recharge,  natural  change  of  the  pathways  of 
circulation  into  the  rock  formations,  losses  of  pressure,  etc.).  If  a 
model  is  considered  to  be  reliable,  it  can  be  used  to  study  how  the 
durability  of  the  resource  changes  depending  on  different  extrac¬ 
tion  rates,  reinjection  temperatures,  wells  siting,  fractures  (also 
induced).  A  productive  strategy  can  be  based  on  the  results  of  a 
model  simulation.  This  is  true  for  both  power  plants  and  thermal 
energy  direct  uses.  Some  general  and  macroscopic  aspects  of  the 
application  of  numerical  simulation  to  the  geothermal  energy 
study  are  listed  here: 

•  equations  describing  the  phenomena  considered  (circulation, 
energy/mass  transport): 

•  estimation  of  the  main  thermophysical  parameters; 

•  boundary  conditions  (BC); 

•  geothermal  potential  assessment: 

•  coupling  between  the  power/thermal  utilization  and  reservoir. 


2.1.  Characterization  of  the  geothermal  source  potential 

Potential  assessment  is  a  fundamental  step  of  a  geothermal 
project  and  its  final  goal  is  the  sustainable  utilization  of  the 
resource.  It  involves  the  complete  characterization  of  the  field, 
energy  stored,  maximum  fluid  rate,  useful  temperatures  and 
chemical  parameters  of  the  fluid  (to  determine  the  minimum 
temperature  for  reinjection).  This  evaluation  is  surely  important 
for  each  kind  of  water  dominant  geothermal  field,  but  mainly  in 
case  of  moderate  temperature  geothermal  fields  in  which  the 
installation  of  an  ORC  plant  is  programmed.  During  the  running  of 
a  plant  the  productivity  (and  wells  deliverability)  can  show  some 
remarkable  variations  (in  terms  of  flow  rate,  fluid  chemistry  and 
specific  enthalpy  of  geofluid).  These  changes  can  be  addressed  to 
reservoir  pressure  decline  because  of  excessive  fluid  production. 
The  reservoir  fluid  pressure  is  also  another  important  parameter, 
linked  both  to  the  productivity  of  the  well  and  the  scaling 
phenomena.  High  pressure  in  the  geofluid  pipes  can  keep  the 
scaling  phenomena  under  controlled  rates. 

2.2.  Reinjection  strategy  and  geofluid  chemistry 

Only  if  reinjection  is  practiced  one  can  say  that  geothermal 
energy  is  being  used  as  a  renewable  energy  source.  The  practice  of 
reinjection  avoids  temperature  and  pressure  decline  in  a  geother¬ 
mal  field.  For  the  binary  power  plants  it  is  a  basic  approach  for 
resource  management  and  it  appears  to  be  compulsory.  The  task  is 
the  optimization  and  enhancement  of  the  durability  of  the 
resource:  the  argument  has  been  diffusely  analyzed  in  [6].  The 
objectives  that  this  practice  has  to  achieve  are  essentially  the 
efficient  restitution  of  the  geofluid  to  the  reservoir  in  order  to 
optimize  the  recharge  in  terms  of  enthalpy  and  flow  rate,  and 
choosing  the  correct  siting  and  depth  of  the  wells  to  guarantee  a 
sustainable  use  of  the  resource.  The  optimum  reinjection  strategy 
is  a  quite  complex  task  and  it  strongly  depends  on  the  type  of 
geothermal  system  (Kaya  et  al.  [9]).  In  general  the  minimum 
temperature  values  are  in  the  range  between  60  °C  and  80  °C. 

The  siting  of  production  and  reinjection  well  and  their  mutual 
interference  are  the  main  issues  of  the  strategy.  To  give  a  trivial 


example:  if  they  are  too  far,  the  recharge  could  occur  in  a  long  time 
interval  (longer  than  the  plant  lifetime).  If  they  are  too  close  a  cold 
fluid  short-circuiting  could  occur.  A  correct  reinjection  strategy  has 
to  take  into  account  the  thermodynamic  and  chemical  problems  of 
the  scaling  phenomena,  which  in  many  cases  increase  with  low¬ 
ering  of  the  temperature.  Early  study  of  the  geothermal  system  is 
fundamental,  in  order  to  avoid  the  worst  consequences  of  fouling, 
corrosion  and  clogging  of  the  parts  of  the  plant,  of  the  pipings,  and 
the  “tapping”  of  the  wells. 


3.  Numerical  simulation  of  geothermal  reservoirs  in  the 
integrated  approach  to  sustainable  design 

Numerical  simulation  of  geothermal  reservoirs  allows  us  to 
understand  the  hydrogeological  behavior  and  heat  transport  into 
the  reservoir  under  a  certain  utilization  rate.  It  is  possible  to  study 
the  geothermal  reservoir  by  solving  the  balance  equations  of  mass, 
momentum  and  energy  in  the  particular  volume  in  which  hydro- 
thermal  circulation  of  fluid  occurs.  The  hydrogeological  issues 
linked  to  the  geothermal  exploitation  (knowledge  of  the  geological 
structures  and  of  the  groundwater  system)  must  be  connected 
with  the  engineering  tasks  of  the  design  and  optimization  of  the 
energy  conversion  system.  The  constitutive  laws  are  peculiar  for 
each  kind  of  reservoir,  while  the  numerical  analysis  issues  are  also 
important  in  order  to  achieve  a  reliable  solution.  The  development 
of  the  numerical  model  itself  has  to  follow  two  main  directions: 

(1)  The  unperturbed  (or  undisturbed)  natural  state. 

(2)  The  utilization  scenarios  (during  the  exploitation). 

Different  phases  in  the  development  of  the  model  can  be 
identified.  A  first  “block-structure”  has  to  be  built  together  with 
the  dataset  of  the  parameters  which  best  fit  what  it  is  expected  by 
the  conceptual  model.  This  first  step  model  should  reproduce: 


Fig.  2.  Conceptual  scheme  for  the  elaboration  of  a  numerical  model  of  a  geother¬ 
mal  reservoir. 
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•  Geological  structure  of  the  reservoir. 

•  Geometrical  features  of  wells  and  fractures. 

•  Hydraulic,  Thermal,  Mechanic  and  Chemical  conditions 

(HTMC). 

3.1.  Conceptual  model  and  numerical  model 

The  conceptual  scheme  for  the  realization  and  simulation  of  a 
model  of  a  geothermal  reservoir  is  represented  in  Fig.  2  (modified 
after  [14]).  If  a  numerical  model  is  realized  in  a  multidisciplinary 
team  or  environment  it  happens  that  who  is  in  charge  of  building 
up  the  model  is  not  totally  aware  of  the  conceptual  model  of  the 
field.  It  then  becomes  necessarily  a  work  of  a  team  in  which 
different  geothermal  backgrounds  and  experiences  are  shared.  An 
important  step  of  “interpretation”  and  appropriate  data  collection 
must  be  pursued  firstly,  for  the  elaboration  of  the  conceptual 
model  (Huenges  [15]).  This  first  step  model  should  reproduce: 
geological  structure  of  the  reservoir;  geometrical  features  of  wells 
and  fractures;  Hydraluic,  Thermal,  Mechanic  and  Chemical  condi¬ 
tions  (HTMC). 

The  model  should  then  pass  a  further  step  of  calibration  and 
refinement.  It  is  an  iterative  process,  in  which  the  parameters  and 
boundary  conditions  should  be  adjusted  according  to  the  con¬ 
ceptual  and  physical  model  previously  elaborated  and  to  the 
uncertainty  level  of  part  of  the  database  used.  Some  thermophy¬ 
sical  properties  (permeability,  thermal  conductivity,  porosity, 
specific  heat  capacity)  change  with  siting,  depth  and  hydrothermal 
alteration.  Their  value  can  be  adjusted  and  the  mesh  can  be  refined 
at  this  step. 

A  first  model  of  the  unperturbed  state  is  the  result  of  the 
attempts  described  here.  It  should  be  run  for  a  long  simulation 
time  interval  (105-106  years)  in  order  to  verify  the  model  stability 
and  convergence.  In  case  of  strong  uncertainty  about  both  the  heat 
transport  phenomena  and  geophysical  data,  simple  2D  models  (or 
lumped  parameters  models)  could  be  firstly  run.  Exploitation  and 
energy  utilization  scenarios  can  be  then  run  starting  from  the 
unperturbed  state  simulation  as  initial  conditions.  The  renew- 
ability  assessment  and  durability  of  the  resource  have  to  be  results 


of  the  scenarios  simulation.  In  particular,  temperature  and  pres¬ 
sure  should  be  kept  stable  in  the  reservoir  as  much  as  possible.  If 
chemical  properties  and  saturation  curve  of  the  specific  geofluid 
mixture  are  known,  scaling  and  chemical  deposition  phenomena 
can  be  also  introduced  in  the  calculations,  in  a  multi-physics 
simulation  environment.  The  models  can  couple  different  trans¬ 
port  equations,  referring  to  mass,  heat  and  chemicals.  It  is 
important  to  remark  that  in  the  results  of  the  simulation  some 
effectively  useful  data  for  the  plant  design  must  be  extracted.  It  is 
evident  that  some  of  the  geophysical  and  general  results  are  not 
directly  necessary  or  relevant  for  who  is  in  charge  of  designing  the 
plant.  A  close  interaction  between  who  elaborates  the  numerical 
model  and  who  designs  the  utilization  plant  would  be  needed 
(according  also  to  the  considerations  and  remarks  about  inter¬ 
disciplinary  work  and  sustainability  assessment  discussed  in  the 
previous  chapters). 

From  a  practical  point  of  view,  building  a  numerical  model  of  a 
geothermal  reservoir  is  not  that  easy.  With  respect  to  other 
engineering  or  science  systems,  a  reservoir  cannot  be  practically 
measured  in  all  its  features.  Its  evolution  can  be  studied,  and  a 
model  can  represent  a  way  of  behavior  prediction  only  if  it  is  based 
on  reliable  data.  In  Fig.  3  a  conceptual  scheme  about  how  the 
development  of  numerical  models  can  help  the  sustainability 
assessment  is  shown.  The  size  of  the  geothermal  area  must  be 
known,  the  geometric  domain  must  be  big  enough  to  comprehend 
every  interesting  mass/energy  transfer  boundary  or  spot,  but  also 
small  enough  to  reach  a  proper  calculation  time.  Present  softwares 
and  mesh  generators  allow  to  set  up  very  complicated  geometries 
(also  polygonal),  but  it  is  a  good  point  to  start  with  simple  mesh 
types  and  domain  shapes  (e.g.  quadrangular,  radial).  A  mesh 
refinement  can  be  useful  only  for  the  area  interested  by  mass/ 
energy  transfers,  like  wells,  natural  recharge,  atmospheric 
geothermal  manifestation.  Anyway  the  problem  of  calculation 
time  saving  must  be  considered,  as  it  can  occur  for  the  operator 
to  think  about  a  massive  refinement,  on  domains  which  are  too 
big.  Once  the  geometry  and  mesh  have  been  set  up,  the  equation 
parameters  (and  constitutive  law)  must  then  be  considered.  The 
constitutive  law  (e.g.  Darcy  Law)  is  typical  for  each  phenomena, 


Fig.  3.  Information  flow  and  prospect  for  the  sustainable  assessment  through  the  use  of  the  numerical  simulation  of  the  reservoir. 
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the  conceptual  model  must  arrange  and  properly  describe  the 
transport  phenomena  occurring  in  the  geological  system.  Thermo¬ 
physical  parameters  should  then  be  put  into  the  simulation.  It  is  a 
big  problem,  as  the  hypothesis  of  homogenous  material  is  useful 
but  it  has  some  practical  limitations.  One  single  value  of  a 
significant  parameter  like  porosity  or  permeability  is  assigned  to 
a  layer,  or  element  (maybe  with  irregular  shape).  Often  one  has  no 
possibility  of  assigning  tensors  or  directional  parameters  (e.g. 
preferential  directions  for  circulation),  as  they  derive  from  very 
accurate  measurements  and  interpretations.  Calibration  is  a  good 
tool,  as  it  allows  to  address  parameter  values  matching  with  real 
temperature  and  pressure  data.  Anyway  it  is  possible  only  when 
reliable  measured  data  are  available.  The  main  parameters  to  be 
assigned  are  essentially: 

•  porosity  (<j6,  defined  as  the  ratio  between  the  voids  and  the  total 
volume  considered)  and  effective  porosity  (taking  into  account  all 
the  “interconnected"  voids,  allowing  then  the  fluid  circulation); 

•  permeability  (k),  being  property  of  both  the  rock  and  the 
circulating  fluid,  giving  an  idea  of  the  productivity  of  such  a 
rock  formation; 

•  density  (p,  its  distribution  derives  from  accurate  measurements 
and  usually  it  is  approximated  with  a  single  value  for  each  rock 
formation); 

•  thermal  conductivity  (of  both  rock  and  fluid), 

•  specific  thermal  capacity  (of  both  rock  and  fluid). 


The  simulation  is  then  run,  and  several  well  known  algorithms 
can  be  used.  The  mathematical/engineering  background  is  also  a 
part  of  this  very  complex  process.  As  the  task  is  the  sustainability 
of  a  project,  a  proper  time  scale  should  be  identified,  considering 
both  economic  lifetime  (20  up  to  50  years)  and  natural  (or 
induced)  renewability  of  the  resource  itself.  Usually  the  steady 
(unperturbed)  state  is  also  simulated,  to  reach  the  conditions 
before  the  exploitation  (in  case  of  new  fields).  In  this  case  the 
time  scale  is  very  large  (104  up  to  107  years).  Large  part  of  the 
models  are  used  to  simulate  industrial  scenarios  of  power/heat 
production  and  study  the  reservoir  response  during  time.  Different 
strategies  can  be  adopted,  in  order  to  keep  the  utilization  feasible 
and  sustainable.  It  is  very  important  to  define  the  size  of  the 
installation  (power  plant,  district  heating)  only  after  the  simula¬ 
tion  and  complete  resource  assessment,  this  being  a  result  of  a 
sustainable  exploitation  strategy,  not  the  starting  point  (Fig.  3). 
Oversizing  or  resource  fast  depletion  (reservoir  cooling,  wrong 
reinjection)  can  be  some  consequences  of  incorrect  sizing. 


3.2.  Mathematical  and  physical  background 


the  Darcy  fluid  velocity  q  is  defined  as  (Saeid  et  al.  [18]) 

k 

q=-(Vp-pgVz)  (1) 

P 

where  k  is  the  intrinsic  permeability,  being  proportional  to  the 
hydraulic  conductivity  I<  according  to  the  definition  (in  [18]) 

K=k£  (2) 

P 


In  which  p  is  the  density,  g  is  the  acceleration  of  gravity  and  p  is  the 
viscosity,  p  is  the  pressure  and  z  is  the  vertical  coordinate.  One 
thing  to  be  considered  is  that  in  case  of  anisotropy  an  hydraulic 
conductivity  tensor  I<  can  be  defined 
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This  element  is  important  in  order  to  consider  the  complexity  of 
the  phenomena  involved.  Average  scalar  values  of  the  hydraulic 
conductivity  or  permeability  are  often  assigned,  but  the  real  behavior 
of  rock  formations  can  be  very  different  from  an  averaged  parameter. 
For  example,  porosity  also  changes  with  pressure  and  depth,  it  is 
important  to  understand  which  assumption  must  be  done  when 
assigning  an  average  value  to  an  extended  rock  body  (having  for 
example  an  unknown  fractures  field).  The  fractures  are  particularly 
difficult  to  reproduce  in  a  numerical  model,  their  accurate  description 
is  important  if  their  spatial  extension  is  to  be  compared  to  the  rock 
formation  size,  otherwise  average  properties  can  be  considered.  In  the 
most  common  reservoir  simulation  softwares,  like  for  example 
TOUGH2  [19],  it  is  also  possible  to  define  different  models  for  the 
relative  permeability  ( R ()  as  a  function  of  the  liquid  saturation  in  the 
mixture  (geothermal  fluid)  (Petrasim  Manual  [20]). 

Let  us  now  consider  a  single  phase  (liquid)  flow  in  a  geothermal 
system.  According  to  Rybach  and  Muffler  [17]  the  following 
assumptions  have  to  be  considered: 

•  rock  matrix  is  homogeneous  and  isotropic,  in  particular  refer¬ 
ring  to  porosity,  permeability  and  thermal  conductivity  (sup¬ 
posed  to  be  independent  by  temperature); 

•  incompressible  fluid;  with  density  (p)  and  kinematic  viscosity 
(u)  dependent  by  temperature  according  to  the  laws: 

p  =  p0[\-a(T-T0)-P(T-T0)2}  (4) 

v=v0a(T)  (5) 

with  po,  v0,  a,  p  and  T0  being  opportune  constants,  while  o(T)  is  a 
function  of  T.  Pressure  work  and  dissipations  due  to  viscosity  can 
be  neglected,  internal  energy  of  liquid  (1)  and  solid  matrix  (r)  being 

E/ =  ci,voi(T —  T0)  (6) 


From  a  mathematical  point  of  view,  the  studies  about  transport  in 
reservoirs  or  aquifers  were  born  in  the  field  of  mining  engineering,  oil 
and  gas  and  applied  thermodynamics.  The  task  is  to  understand  the 
response  of  a  porous  or  fractured  system  when  a  hydraulic  gradient  is 
applied.  Anyway  the  application  of  some  experimental  concepts  to  the 
geothermal  fields  is  not  a  trivial  task.  The  hypothesis  and  the  real 
situation  of  the  field  should  be  compared  very  carefully. 

The  basic  calculations  executed  by  the  softwares  for  numerical 
simulation  substantially  deal  with  the  resolution  of  the  flow  into 
porous  and/or  fractured  media.  Equations  of  conservation  of  the 
following  properties  are  solved:  mass,  momentum,  energy,  con¬ 
centration  of  pollutant  (or  chemicals)  dissolved.  Some  constitutive 
laws  must  also  be  implemented,  dealing  with  the  particular 
phenomenon  involved.  Geothermal  reservoirs  fluid  flow  can  be 
generally  studied  according  to  the  Darcy  Law  about  porous  media, 


Er  =  crvo,(T-T0)  (7) 

with  ctvo(  and  cr-vo;  specific  volumetric  heat  capacities  of  rock  and 
liquid.  Under  these  assumptions  the  balance  equations  of  mass, 
momentum  and  energy  can  be  written  respectively  as  (according 
to  [17]): 


Vq  =  0 


YP)+a(T-T0) 

Po) 


1+qr-To) 

a 


t+v-f-t  =  o 


(8) 

(9) 


,dT 


[(1  —  Pl/pCy  +  i/ipoCv  —  +P()CV  q  '  VT  —  V  ■  (kmVT) 


(10) 


km  being  the  thermal  conductivity  of  the  mixture  solid  liquid. 
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A  double  phase  system,  in  which  the  vapor  phase  is  also 
considered,  is  described  by  similar  equations  [17],  The  following 
assumptions  about  porosity  have  to  be  done  in  this  case: 

•  porosity  0  only  depends  on  local  pressure  of  fluid; 

•  liquid  and  vapor  phases  are  in  local  thermal  equilibrium. 


The  equation  of  conservation  (respectively)  of  mass,  momen¬ 
tum  (liquid  and  vapor)  and  energy  can  be  then  written  as  follows: 


^(<pp)  +  V  ■  (Pi~q  ,  +Pv~qv)  =  0 

HD 

— A/  K  /  — >-\ 

q  i  =  — —  ■  [yp-pig) 

(12) 

— >•  Ay  K  — y 

?»= - '(Vp-pvg) 

(13) 

—  [(1  —  </>)prEr  +  tppE]  +  V(/jj£|  q  i+pvEy  q  v+p  q  i+p  q  v) 

=  V(fcmVT) ■+•(/>(  q  i~EPv  q  v) '  if  (14) 

where  E  is  the  internal  energy  of  the  liquid-vapor  mixture,  k  is 
the  permeability  tensor  and  is  the  relative  permeability  of  the 
i-th  phase  ( i=l ,  g  liquid  or  vapor).  Some  fundamental  differences 
between  the  flow  in  porous  and  coherent  media  than  in  fractured 
media  must  then  be  taken  into  account  ([7,17]): 

(a)  permeability  induced  by  fracturing  is  generally  higher  than 
average  formation  permeability: 

(b)  fracturing  permeability  is  usually  anisotropic  (it  depends  on 
fracturing  preferential  direction): 

(c)  the  permeability  due  to  fracturing  is  considerably  more 
dependent  on  pressure  and  tension  field  in  the  rock  with 
respect  to  rock  matrix  permeability. 


In  the  recent  paper  by  Zeng  et  al.  [21]  a  discussion  about  the 
methods  for  fractured  systems  simulation  is  given.  Geothermal 
fracture  systems  can  be  represented  essentially  in  two  ways: 

(1)  discrete  fracture  network  method  (fracture  orientation,  spa¬ 
cing  and  mechanical  properties): 

(2)  equivalent  continuous  porous  media  method  (the  fracture  system 
is  seen  as  an  equivalent  continuous  media,  with  similar  meth¬ 
odologies  as  for  double  porosity  and  double  permeability). 

In  Fox  et  al.  [12]  an  analytical  and  numerical  model  for 
fractured  reservoirs  (multi-fractured  EGS)  is  presented  and  simu¬ 
lated.  In  the  following  Fig.  4  an  example  of  simplified  scheme  of 
fluid  circulation  into  the  fractures  is  shown  ( b  being  the  fracture 
size  and  xs  the  fractures  mutual  spacing).  There  is  also  the 
possibility  of  coupling  thermal-flow  problems  with  chemical  or 
geomechanics  problems  into  the  same  numerical  simulation.  It  is 
the  case  of  some  models  presented  in  literature,  like  for  example 
the  case  of  a  fully-coupled  flow  and  geomechanical  model  by  Hu 
et  al.  [22], 

3.3.  Boundary  and  initial  conditions 

The  structural-physical  model  has  to  be  set  with  the  boundary 
and  initial  conditions,  and  also  with  constraints  that  can  be  useful 
for  the  results  stability.  The  thermal  boundary  conditions  (BC) 
usually  represent  the  heat  geothermal  source  entering  the  reser¬ 
voir:  heat  flow  from  the  bottom,  fixed  temperature  values  at 
bottom/top  or  intermediate  layers  and  adiabatic/impermeable 
conditions.  BCs  usually  also  represent  the  natural  manifestations, 
natural  recharge,  lateral  or  regional  flows,  wells  withdrawal/ 
reinjection  in  the  aquifers  and  hydraulic  head  both  for  the 
hydrological  problem  and  for  the  heat  transfer.  The  BC  kind  are 
similar  when  considering  flow  or  heat  transport.  Fig.  5  provides 


Fracture 


Fig.  4.  Conceptual  scheme  for  multi-fracture  vertical  geothermal  system  (xs  is  the  spacing  between  fractures,  b  is  the  fracture  aperture). 
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Fig.  5.  Example  of  boundary  conditions  (both  thermal  and  mass  conditions)  in  a  3D 
numerical  model  grid. 

example  of  boundary  conditions  in  a  3D  grid.  Mainly  four  types  of 
BC  are  used: 

( 1 )  First  kind  (or  Dirichlet )  condition  -  along  a  border  or  a  boundary 
the  hydraulic  head  or  temperature  are  assigned. 

(2)  Second  kind  (or  Neumann )  condition  -  along  a  border  or  a 
boundary  the  fluid  or  heat  flux  is  assigned. 

(3)  Third  kind  (or  Cauchy )  condition  -  transfer  coefficients  are  used 
particularly  for  the  hydraulic  head. 

(4)  Fourth  kind  condition  -  single  well  or  singular  point  source, 
typically  used  for  wells  conditions  (extraction  or  reinjection, 
according  to  the  sign  convention),  implementing  Dirac  S 
function. 

In  particularly  dry  reservoirs,  conduction  is  the  prevalent 
mechanism  of  heat  transfer.  In  hydrothermal  aquifers  and  tradi¬ 
tional  geothermal  systems  convection  flow  also  contributes  to  the 
mass/heat  transport  phenomena.  Thermophysical  parameters 
database  are  also  available  in  the  most  used  codes.  Anyway,  as  it 
is  stated  in  this  work,  a  characterization  of  the  parameters  should 
by  site-dependent  in  order  to  obtain  reliable  and  physically 
consistent  results. 

The  initial  conditions  are  usually  the  thermal  gradient,  pressure 
distribution  in  the  domain  and  hydraulic  head  levels  of  rivers  or 
reservoirs.  To  hold  the  results  range  near  a  specific  value,  con¬ 
straints  about  max/min  fluid  rate,  pressure  or  heat  flow  can  be  set. 
Fig.  6  summarizes  the  various  steps  described  here:  first  the 
definition  of  geological  elements  and  dimension  of  the  reservoir, 
then  spatial  discretization  and  materials  calibration,  finally  the 
definition  of  boundary  conditions,  initial  conditions  and  definition 
of  temporal  domain. 

3.4.  Numerical  integration  of  the  equation  systems 

The  integration  of  all  the  interdisciplinary  inputs  and  proce¬ 
dures  is  the  most  challenging  and  crucial  part  of  a  modeling 
process.  An  important  issue  of  this  process  deals  with  the  quality 
of  information  and  data  flowing  through  the  starting  phase  and 
the  simulation  itself.  Numerical  simulation  must  be  treated  and 
used  as  an  iterative  process,  continuously  changing  and  improving, 


Geology  -  Spatial  extension  -  Rock  materials 


* 


3D  features  -  Spatial  discretisation  - 
Materials  calibration 


Fig.  6.  Overall  graphic  view  of  the  numerical  simulation  process  (the  second 
vertical  geological  section  is  from  [23]). 


as  the  information  flow  goes  both  ways  (Ungemach  et  al.  [14]). 
A  good  and  reliable  model  is  often  an  object  of  discussion  and 
reinterpretation.  There  are  different  techniques  for  space  discreti¬ 
zation:  finite  difference,  finite  volume  methods,  as  well  as  finite 
element  methods.  Different  numerical  integration  techniques  are 
implemented  in  codes  and  softwares.  TOUGH2  [19]  for  example, 
uses  a  finite  difference  (time  fully  implicit)  discretization,  while  in 
Feflow  [24]  finite  element  techniques  are  used.  Other  softwares 
are  used  in  literature  and  in  industry.  Many  softwares  have  been 
developed  and  used  by  Research  Institutes  or  Universities.  In  those 
softwares  a  lot  of  the  most  known  resolution  algorithms  (from 
numerical  analysis  and  calculus)  are  implemented.  The  mesh 
refinement  is  a  fundamental  instrument  that  can  be  adopted  to 
improve  the  analysis  and  optimize  the  computational  tasks  (con¬ 
centrating  for  example  the  mesh  number  in  the  wells  area  or  along 
the  faults).  Different  techniques  can  be  adopted  for  modeling  of 
the  faults,  but  the  accuracy  about  the  data  in  input  (upflow / 
downflow,  fluid  rate,  permeability,  thermal  anomaly)  for  these 
structures  has  to  be  very  high  to  achieve  reliable  results. 

All  the  most  used  softwares  are  multipurpose,  involving  the 
possibility  of  simulating  different  types  of  diffusion  phenomena. 
Pre  and  post-processors  are  typically  embedded  in  commercial 
softwares,  so  that  graphical  interface  and  elaboration  of  the  data 
can  be  easily  carried  out.  The  various  design  tools  (also  for  shallow 
systems)  are  mainly  based  on  the  line-source  model  and  Eskilson's 
approach,  and  the  duct  ground  storage  model  DST  (Haberl  and  Do 
[25]).  When  coupling  of  more  phenomena  is  implemented  it  can 
give  a  very  useful  contribution  for  the  scaling  phenomena  limita¬ 
tion  (concentration,  inflow-outflow,  and  variation  during  with¬ 
drawal-reinjection  must  be  known). 
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Fig.  7.  Typical  dotted  breakthrough  temperature  profile:  (a)  single  extraction  well;  (b)  and  (c)  double  extraction  well. 


3.5.  Limitations  and  criticalities  of  numerical  simulation  of 
geothermal  reservoirs 

Notwithstanding  its  strategic  importance  it  is  clear  that  numer¬ 
ical  simulation  of  geothermal  reservoirs  has  some  important 
limitations  (which  are  treated  and  discussed  in  this  work): 

•  it  strongly  depends  on  the  reliability  and  accuracy  of  the  data; 

•  numerically  stable  models  can  be  physically  inconsistent. 

The  first  limit  can  be  also  expressed  by  the  principle  “trash  in¬ 
trash  out”.  It  must  be  clear  which  is  the  physical-numerical 
problem  to  recourse  numerical  simulation  for.  One  should  evaluate 
if  the  numerical  simulation  is  the  more  appropriate  instrument  to 
face  a  specific  problem.  Usually  a  problem  can  be  simplified  in  a 
proper  way  to  be  solved  according  to  calculus  or  numerical 
analysis  without  using  dedicated  softwares  and  complicated 


geometric  domains.  Reservoir  model  simulation  has  to  be  pursued 
only  if  it  is  the  most  suitable  and  appropriate  way  to  elaborate  a 
design  strategy.  “Lumped  parameters”  models  can  be  very  useful 
for  some  medium-low  temperature  fields,  considering  plain  litho¬ 
logical  layers.  Sometimes,  particularly  for  linear  and  simple  pro¬ 
blems  they  can  be  satisfactory,  in  spite  of  more  sophisticated 
elaborations. 

One  possible  risk  is  to  start  “asking  too  much”  or  “asking  too 
bad  response”  to  the  numerical  models,  making  them  “too  much” 
or  “too  less  powerful”.  For  example,  starting  from  the  same 
geological  features  of  a  field,  a  model  can  give  different  results 
depending  on  the  resolution  of  the  spatial  distribution  of  the  data. 
The  numerical  analysis  could  be  important  in  order  to  define  the 
exploitation  scenario  (Figs.  7  and  8)  or  the  reinjection  strategy. 

Some  typical  problems  due  to  incorrect  initial  characterization 
of  the  resource  can  also  be  faced  with  an  appropriate  model 
simulation,  they  can  deal  with:  oversizing  of  the  plant  (leading  to 


A.  Franco,  M.  Vaccaro  /  Renewable  and  Sustainable  Energy  Reviews  30  (2014)  987-1002 


995 


excessive  extraction),  scaling  (causing  corrosion,  productivity 
drop,  diameter  reduction)  and  wrong  reinjection  strategy  (losses 
of  fluid  or  cooling  of  the  reservoir). 

4.  Numerical  simulation  of  geothermal  reservoirs:  a  review  of 
case  studies  referred  to  experimentally  investigated 
geothermal  fields 

In  this  part  of  the  paper  various  examples  of  numerical  simula¬ 
tions  of  geothermal  reservoirs  are  analyzed  and  critically  reviewed. 
The  point  of  connection  among  the  various  reviewed  cases  is  the  fact 
that  they  are  referred  to  very  well  known  geothermal  fields  for  which 
a  lot  of  experimental  data  are  also  available.  Each  simulation  has  its 
own  peculiarity,  concerning  with  dimensional  scale,  enthalpy  level  of 
the  reservoir,  fluid  rate  extracted/reinjected,  and  software  used.  A 
summary  of  the  meaningful  data  is  also  given.  All  the  cases  described 
are  different  for  various  reasons.  First  of  all  for  the  typology  of 
geothermal  field:  from  medium  enthalpy  water-dominant  field  to 
dry-steam  dominant  field.  The  differences  between  the  models  deal 
with  simulation  domains  (size  ranges  from  some  1cm  up  to  about 
100  km),  scenarios  simulated  (unperturbed  or  exploitation)  and 
software  used. 

One  concept  has  to  be  emphasized  in  this  review:  the  strong 
dependence  of  the  results  of  the  numerical  analysis  on  the  quality 
of  the  inputs  and  the  difficulty  that  would  be  afforded  in  realizing 
the  models.  First  of  all,  the  data  and  the  geo-thermo-physical 
parameters  necessary  are  not  always  available  or  measurable. 

In  the  next  two  sections  some  cases  will  be  analyzed:  the  review  is 
organized  as  follows:  in  the  first  part  are  a  group  of  cases  related  to 
geothermal  field  for  which  a  lot  of  experimental  data,  boundary 
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Fig.  8.  Graphical  results  about  a  simplified  single  well  extraction  system  (also  fluid 
velocity  vectors  and  iso-temperature  vertical  section  are  shown). 


conditions  are  available  and  it  is  possible  to  reproduce  the  results  of 
the  simulation  (Table  1).  Each  case  is  discussed  in  detail  in  a  single 
subsection.  A  second  group  concerns  cases  available  in  the  literature 
and  refers  to  less  investigated  (and  simulated)  geothermal  systems. 
The  different  techniques  of  numerical  resolution  of  the  equations  in 
the  model  are  also  considered  here. 

4.1.  Momotombo  reservoir  (Nicaragua) 

The  first  case  analyzed  is  the  geothermal  field  of  Momotombo 
in  Nicaragua  for  which  a  meaningful  set  of  data  is  available  from 
literature.  This  case  study  has  been  chosen  to  remark  how 
important  the  characterization  and  assessment  of  the  geothermal 
potential  is  to  undertake  a  correct  industrial  exploitation  of  a 
reservoir.  This  analysis  is  based  on  the  literature  data  available, 
mainly  from  the  work  of  Porras  et  al.  [26-28], 

Overestimation  of  the  geothermal  potential  (and  progressive 
plant  oversizing)  brought  to  a  gradual  and  progressive  impover¬ 
ishment  of  the  resource  in  terms  of  temperatures,  pressures  and 
geothermal  fluid  rates.  Since  the  initial  years  of  industrial  interest 
(1983-1989)  the  production  decreased  and  different  problems 
arised.  A  three-dimensional,  porous-media,  numerical  model  of 
the  system  was  developed  and  calibrated  and  it  was  utilized  to 
study  the  response  of  the  geothermal  reservoir  under  different 
fluid  production  and  injection  scenarios.  The  initial  hypotheses 
are:  total  reinjection  of  the  fluid  extracted;  reinjection  at  100  °C 
and  pressure  5  bars;  time-constant  fluid  rates. 

The  domain  considered  in  this  case  is  3.1  x  2.4  km2  with  a  total 
depth  of  -3000  m,  divided  into  9  layers.  The  whole  domain  is 
built  with  972  blocks  (minimum  dimensions  200  x  200  m2,  max¬ 
imum  dimensions  600  x  700  m2).  The  total  number  of  materials 
considered  is  18.  The  most  permeable  layer  is  5  x  lCrt11  m2,  in  a 
depth  interval  between  -150  m  and  -450  m.  Four  different 
scenarios  of  exploitation,  dealing  with  different  extraction/reinjec¬ 
tion  strategies  are  proposed.  In  [26-28]  the  authors  state  that 
satisfactory  agreements  were  obtained  between  measured  and 
computed  discharge  enthalpies  and  flow  rates,  for  most  of  the 
shallow  wells.  The  model  also  qualitatively  reproduced  the  pres¬ 
sure  drawdown  measured  in  selected  wells. 

4.2.  Ngatamariki  geoihermal  field  (New  Zealand) 

The  Ngatamariki  geothermal  field  is  located  17  km  from  Taupo 
(North  Island,  New  Zealand).  It  is  one  of  the  several  high  enthalpy 
geothermal  systems  within  the  Taupo  Volcanic  Zone  (they  are 
more  than  20).  Exploration  wells  were  first  drilled  in  1985-1986, 
and  Mighty  River  Power  then  drilled  further  wells  in  2008-09.  The 


Table  1 

Summary  about  the  various  numerical  model  analyzed:  data  referred  to  experimentally  investigated  reservoir. 


Reservoir 

Software 

Domain  size  and  blocks  no. 

Mean  res.  T  [  C] 

Fluid  flow  rate 

Momotombo 

(Porras  et  al.,  [26-28]) 

TOUGH2 

3.1  x  2.4  km2;  3  km  depth  ;  972 
blocks 

240-340 

—  357  kg/s;  13  production  wells; 

8  reinjection  wells  (power  units:  32  MW) 

Ngatamariki 

(Burnell  and  Kissling  [29]) 

TOUGH2 

10.5  x  11  km2;  —5  km  depth;  24128 
blocks 

80  (scenario  A);  120 
(scenario  B) 

—  695  kg/s  (2  scenarios) 

Larderello  1  (Romagnoli  et  al.  [30], 
Barelli  et  al.  [31],  Arias  et  al.  [32]) 

TOUGH2 

70  x  70  km2;  —7.5  km  depth; 

-10000  blocks 

200-300  C 

1300  kg/s  (average  historical  data) 

Larderello  2  (Della  Vedova  [33]) 

SHEMAT 

42  x  26  km2;  10  km  (total  thickness) 

200-300  C 

(Only  natural  state  simulation) 

Wairakei  (O'Sullivan  et  al.  [34,35]; 
Mannington  et  al.  [36,37]) 

Non-commercial 
codes;  TOUGH2 

30  x  30  km2;  3.4  km  (total  thickness); 
8055  blocks  (2006  model) 

250  °C  (Wairakei-Te  Mihi) 

— 1460  kg/s  (historical  average);  >  200 
wells  (57  years)  (power  units:  >  200  MW) 

GroR  Schonebeck 

(Blocher  et  al.  [39]) 

FEFLOW 

—4.8  x  5.5  km2;  —0.6  km  depth; 
489591  prism  elements;  254744 
nodes 

150  °C  (125  °C  after 

30  years  simulation) 

—21  kg/s  doublet  of  wells  from  the  same 
drilling  area 

Mt.  Amiata  (Barelli  et  al.  [40]) 

TOUGH2 

1100  km2;  >  5.7  km  (total  thickness) 

150-220  °C  (first  shallow 
res.);  300-350  °C(deep  res.) 

(Natural  state  simulation;  interaction 
between  reservoirs) 

Dubti  (Tendaho  Rift) 

(Battistelli  et  al.  [41,42]) 

TOUGH2  / 

EWASG 

—2.5x3  km2;  >2  km  depth 

245-290  °C 

140  kg/s 

996 
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company  has  planned  the  installation  of  overall  130  MW  electric 
power  output,  from  a  condensing  unit  and  an  ORC  unit  integrated 
in  a  combined  cycle  configuration  (to  be  first  run  in  2013).  The 
data  about  this  geothermal  field  are  available  from  a  group  of 
industrial  reports  by  Burnell  and  Kissling  (2009),  in  particular  the 
one  about  the  reservoir  model  [29],  in  which  a  very  complete 
conceptual  model  can  be  found.  The  geothermal  reservoir  is 
composed  of  3  parts:  a  shallow  aquifer  (50-100  m  deep);  an 
aquitard  in  the  rock  formation  of  Huka  Falls  (50-400  m  deep), 
between  the  surface  aquifer  and  the  intermediate  one;  and  an 
intermediate  aquifer,  between  Huka  Falls  rock  formation  and  the 
clay  cap.  The  grid  has  globally  26  horizontal  layers  (of  different 
thickness),  subdivided  into  928  elements,  reaching  a  total  number 
of  24128  blocks.  The  surface  extension  of  the  model  is  132  km2.  An 
inlet  of  fluid  flow  (70  kg/s)  is  considered,  entering  from  the  bottom 
of  the  model,  at  a  specific  enthalpy  of  1450  kj/kg,  with  an  heat 
flow  of  101.5  MW.  From  the  basement  a  constant  heat  flow  of 
11  MW  is  considered.  As  constraints,  areas  with  known  flow 
direction  and  impermeable  formations  are  simulated  using  fixed 
state  (Tgeo,  Pgeo)  conditions.  The  validation  of  the  model  is  based  on 
temperatures  and  pressure  data  measured  in  different  wells  and 
on  the  calibration  of  other  parameters  like  porosity,  permeability, 
upflow  mass  and  enthalpy,  hydraulic  connections  between  the 
deep  reservoir  and  the  groundwater  system. 

4.3.  Larderello  geothermal  field  (Italy),  2010  model 

Larderello  field  (Italy)  is  one  of  the  most  anciently  known  and 
studied  geothermal  areas  of  the  world.  This  field  has  been  widely 
drilled  and  developed,  with  an  almost  100-year-old  history  in 
geothermal  energy  utilization  for  power  purposes.  Average  fluid 
production  in  the  Larderello  field,  after  the  most  recent  explora¬ 
tions  and  improvements  is  now  about  3700  t/h.  The  exploration 
extended  in  the  early  1950s  to  the  near  Travale  field  (10-15  km  SE 
of  Larderello),  which  has  now  increased  its  fluid  production  to  an 
average  value  of  1000  t/h.  For  this  reason,  when  talking  about 
large  scale  model  of  this  geothermal  system,  usually  one  can  talk 
about  Larderello-Travale  geothermal  field. 

A  numerical  model  about  the  field  has  been  recently  realized 
and  improved,  increasing  the  dimensions  of  the  geological  domain 
considered,  by  Romagnoli  et  al.  [30],  Barelli  et  al.  [31],  Arias  et  al. 
[32],  The  model  domain  extent  is  4900  km2  (70  km  each  side), 
with  a  total  thickness  of  7.5  km.  The  grid  is  made  of  10000  cells 
and  16  vertical  layers.  The  geological  scheme  refers  to  five  main 
rock  formations:  clayey-shaley  caprock  (0-500  m),  fractured  car¬ 
bonate  reservoir  (500-1000  m),  metamorphic  reservoir  (1000- 
5000  m)  and  granitic  intrusion  as  heat  source  of  the  system. 
Sixteen  rock  materials  are  considered  and  an  impermeability 
condition  along  the  boundaries  is  assigned.  Fixed  state  (time 
invariant)  conditions  of  temperature  are  considered  at  the  top 
(15  °C,  atmospheric  pressure)  and  at  the  bottom  of  the  producing 
layer  (350-400  °C).  Natural  manifestations  and  cold  inflow  from 
the  shallow  aquifers  are  the  only  interactions  with  the  external 
environment.  The  simulation  of  natural  state  has  been  carried  out 
(millions  of  years  as  temporal  scale)  and  a  simulation  of  the 
exploitation  history  of  the  field  has  also  been  modeled.  The 
historical  data  of  700  wells  have  been  represented  with  20 
“virtual”  wells.  The  conclusion  of  Romagnoli  et  al.  [30]  are  that 
only  few  changes  in  the  conditions  of  the  natural  system  have 
been  caused  by  industrial  development  of  the  area. 

4.4.  Larderello  geothermal  field  (Italy),  2008  model 

A  different  model  of  this  field  has  been  proposed  by  Della 
Vedova  et  al.  [33].  This  model  deals  only  with  the  natural  state  of 
the  geothermal  system,  without  considering  the  industrial 


exploitation.  A  very  large  temporal  scale  is  considered  (8-12 
millions  of  years).  The  extent  of  the  considered  domain  are 
42  x  26  km2,  with  a  total  thickness  of  10  km.  The  total  depth  of 
the  model  is  very  high,  to  include  the  K-horizons  and  the  data 
from  fluid  inclusions.  K-horizons  are  considered  to  be  the  main 
reservoir  bottom,  corresponding  to  the  400  °C  isotherm;  collo¬ 
cated  between  3000  m  in  the  Larderello  zone  and  104  m  in  the 
Travale  zone.  The  numerical  simulator  SHEMAT  has  been  tested 
and  validated  with  the  geophysical  data  available  for  the  upper 
crust.  The  mesh  cell  size  is  1  x  1  x  0.3  km3.  The  upper  surface 
boundary  conditions  are  20  °C  and  0.1  MPa,  impermeable  and 
adiabatic  conditions  are  assigned  at  the  lateral  boundaries.  The 
bottom  boundary  is  assumed  to  be  impermeable  and  at  a  fixed 
temperature  of  400-600  °C.  A  sensitivity  analysis  about  thermal 
parameters  is  also  considered  in  [33],  The  authors  remark  that  a 
lot  of  simulation  have  been  run  to  match  a  composite  target 
function,  due  to  the  uncertainty  about  several  input  data  (geome¬ 
try,  rock  data).  The  work  considered  is  an  example  of  deep  field 
simulation,  oriented  to  the  comprehension  of  the  deep  field 
phenomena  more  than  to  a  sustainable  exploitation  approach. 

4.5.  Wairakei  geothermal  field  (New  Zealand) 

Wairakei  is  another  well-known  geothermal  field  located  in  the 
Taupo  Volcanic  zone  (New  Zealand).  The  extent  of  the  area  is 
almost  25  km2.  A  complete  report  about  the  industrial  history  and 
the  evolution  of  the  numerical  models  can  be  found  in  O'Sullivan 
et  al.  [34]  and  Bixley  et  al.  [35],  Fifty  years  of  activity  on  this  field 
have  been  reached  in  2008.  In  the  period  1958-1990  the  power 
installed  was  increased  up  to  140  MW.  After  1990  the  power 
installed  has  approached  200  MW  (with  the  plant  of  Poihipi, 
55  MW).  The  trend  of  the  specific  enthalpy  of  the  geofluid  in  the 
field  fits  the  evolution  during  the  years  [35]:  rise  of  pressure  due 
to  reinjection  activity;  increase  of  temperature  in  the  wells  in  the 
“Eastern  Borefield”  and  seepage  of  cold  water  in  the  reservoir. 

Several  models  have  been  proposed  and  tested  for  Wairakei 
reservoir  [34],  The  models  have  been  used  to  simulate  different 
scenarios  of  field  development.  The  common  aspects  and  main 
characteristics  of  the  conceptual  schemes  are:  big  fractures  and  faults 
(NE-SW),  increasing  permeability;  two  big  upflow  areas  (260  C  at 
Wairakei,  300  °C  at  Tauhara);  heat  flux  fixed  values  as  BC,  in  the  range 
300-600  MW;  two  regions  at  different  pressures  (high  pressure  at  Te 
Mihi,  low  pressure  at  Southern  Wairakei);  natural  recharge  and 
rainwater  infiltrations  (1000  mm/year,  5%  infiltration).  A  discussion 
about  the  calibration  of  the  model  is  also  reported  in  Mannington  et  al. 
[36,37],  The  code  iTOUGH2  [38],  a  TOUGH2  family  software  for  inverse 
simulation,  has  been  used  to  support  the  calibration  process  of  the 
parameters  and  improve  the  match  to  the  field  data.  Good  matching 
results  have  been  reached  when  the  Tauhara  field  (strictly  connected 
with  Wairakei)  has  been  introduced  in  the  domain  of  the  model.  The 
perspective  is  to  improve  a  27886-blocks  model  (comprehending 
Wairakei-Tauhara-Rotokawa). 

4.6.  Crofi  Schonebeck  reservoir  (Germany) 

A  case  study  about  a  simulation  of  an  Enhanced  Geothermal 
System  (EGS)  in  a  deep  geothermal  reservoir  with  hydraulic 
stimulation  is  discussed  in  Blocher  et  al.  [39].  The  geothermal 
research  site  of  Grol?  Schonebeck  is  located  40  km  north  of  Berlin 
(Germany).  This  case  is  presented  as  an  example  of  well-based 
hydrothermal  problem  with  a  good  match  with  the  sustainable 
energy  exploitation  issue.  A  very  good  approach  both  to  the 
numerical  simulation  and  to  the  accuracy  of  the  data  is  achieved. 

The  software  used  for  the  simulations  is  FEFLOW  [23],  The 
simulation  deals  with  a  doublet  of  wells.  The  reservoir  is  located 
between  -3815  and  -4247  m  below  sea  level,  in  the  Lower 
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Permian  of  the  Northeast  German  Basin.  The  hydraulic  stimulation 
is  a  technology  used  to  improve  the  productivity  of  a  reservoir  by 
inducing  artificial  fractures  through  high  pressure  fluid  injection 
(water,  gel-proppant).  Six  geological  formations  are  considered, 
with  different  lithologies.  Two  sandstones  formations  (between 
-4000  m  and  -4100  m  depth)  are  the  most  appropriate  for 
geothermal  power  production.  Natural  fracture  system  is  studied 
to  analyze  the  conceptual  scheme  of  circulation  of  the  water,  and 
to  develop  the  induced  fractures  layout.  The  wells  have  a  distance 
of  28  m  at  the  surface,  the  injection  well  is  vertical,  while  the 
production  well  is  deviated  to  guarantee  a  sufficient  distance  of 
500  m  between  them  within  the  reservoir. 

An  interesting  aspect  remarked  in  [39]  is  the  dependence  of  the 
hydraulic  properties  with  temperature  and  pressure  in  the  reser¬ 
voir  (density,  permeability,  porosity,  thermal  conductivity)  to  be 
considered  as  part  of  the  model.  Also  the  total  dissolved  solids  and 
chemical  composition  are  considered  in  the  whole  model.  The 
model  has  an  extension  of  4.809  x  5.448  km2,  with  a  depth  of 
almost  0.6  km.  The  grid  is  made  by  489591  prismatic  elements  and 
254744  nodes,  discretizing  27  spatial  layers.  The  induced  fractures 
are  represented  by  2D  quadrilateral  vertical  elements.  The  produc¬ 
tion  well  has  a  fluid  rate  of  75  m3/h  (~21  kg/s)  at  150  °C  with  a 
concentration  of  solids  of  265  g/1.  The  injection  well  has  the  same 
fluid  rate  and  concentration  of  solids  but  the  temperature  is  70  °C. 
The  quasi-stationary  conditions  are  reached  after  almost  4  x  104 
years  of  simulation  time  (hydraulic  head  levels  are  matched).  A  30 
years  simulation  of  the  exploitation  with  the  above  mentioned 
wells  conditions  gives  a  production  temperature  drop  from  150  °C 
to  125.8  °C. 

4.7.  Mt.  Amiata  geothermal  system  (Italy) 

In  a  recent  paper  of  Barelli  et  al.  (2010)  the  Tuscan  geothermal 
system  of  Mount  Amiata  is  considered  [40],  The  fields  of  Bagnore 
and  Piancastagnaio  have  been  explored  and  drilled  (more  than  100 
wells)  for  about  50  years.  In  this  field  two  main  reservoirs  are 
present:  the  first  one  is  in  the  carbonatic  formations,  between  500 
and  1000  m  deep,  at  average  temperature  of  150-220  °C;  the 
second  reservoir  is  in  the  Paleozoic  metamorphic  basement  at 
depths  of  2500-4000  m,  at  temperatures  of  300-350  C.  The  Mt. 
Amiata  Volcanic  Complex  has  a  total  area  of  80  km2.  The  peculiar¬ 
ity  of  the  model  is  not  only  to  analyze  the  exploitation  of  the 
reservoir  but  also  to  verify  the  possibilities  of  interaction  between 
a  phreatic  aquifer  (separated  from  the  shallow  aquifer  by  few 
hundred  meters  of  impermeable  formations)  and  the  geothermal 
system.  Moreover  another  particular  aspect  is  that  the  two  main 
reservoirs  are  characterized  by  gaseous  caps  (gas,  vapor),  in 
structures  named  “traps”.  This  is  a  peculiarity  of  Mt.  Amiata  field, 
occurring  in  the  definition  of  the  pressure  distribution  and  fluid 
circulation  model.  The  numerical  model  considered  has  been 
simulated  with  the  software  TOUGH2  [19],  The  surface  extension 
of  the  model  is  more  than  1100  km2,  with  a  total  thickness 
between  -4  km  and  1.738  km  (a.s.I.,  Mt.  Amiata  top  elevation). 
A  time-constant  heat  flow  (average  400  mW/m2,  with  peaks  of 
600-700  mW/m2)  is  the  bottom  boundary  condition.  The  model  is 
globally  closed,  referring  to  inflow-outflow  of  water.  As  remarked 
in  [40],  the  outputs  match  with  the  field  data. 

4.8.  Dubti  geothermal  field  (Tendaho  Rift,  Ethiopia) 

The  model  of  the  Dubti  geothermal  field  review  is  based  on  a 
paper  by  Battistelli  et  al.  [41  ].  The  Tendaho  Rift  was  identified  as  a 
promising  geothermal  area  since  an  exploration  project  of  the  late 
1960s,  and  early  1970s.  More  recently  (1990s)  a  drilling  and  wells- 
testing  activity  was  carried  out  to  explore  and  assess  the  geother¬ 
mal  resource  present  in  the  area  (central  part  of  Northern  Tendaho 


Rift).  The  drilling  in  the  area  of  Dubti  plantation  confirmed  the 
existence  of  a  shallow  geothermal  reservoir.  The  temperature 
recorded  in  the  drilled  wells  is  in  the  interval  245-270  °C,  while 
the  temperature  of  the  geofluid  rising  in  a  fault  in  the  area  is  about 
290  °C  (natural  manifestations  and  fumaroles  are  present  in  the 
whole  area).  Lots  of  production/injection  tests  have  been  carried 
out  in  the  area,  and  the  paper  cited  is  very  detailed  about  these 
data  and  permeability  distribution.  The  Dubti  fault  is  considered  to 
be  ascending  and  upflow  hypothesis  is  at  the  base  of  the 
conceptual  model.  Horizontal  circulation  (with  eventually  two- 
phase  conditions),  when  crossing  permeable  layers,  is  also  con¬ 
sidered.  For  the  numerical  model  the  software  TOUGH2  [19]  has 
been  used,  implementing  the  EWASG  equation-of-state  module 
developed  for  water,  sodium  chloride  and  C02  mixtures  (Battistelli 
et  al.  [42]).  The  simple  3D  model  is  extended  about  7.5  km2  with  a 
total  thickness  of  about  2  km.  Natural  meteoric  recharges  are 
considered  (coming  from  the  Ethiopian  Plateau)  together  with 
horizontal  and  sub-vertical  flows.  The  assumed  initial  temperature 
of  the  hot  upflow  is  290  °C.  Model  results  and  further  drillings 
confirmed  the  presence  of  a  hot  fluid  circulation  zone  at  depths 
between  250  and  500  m  in  the  Dubti  area,  also  confirming  the 
existence  of  more  permeable  zones  along  the  Dubti  fault.  A 
possible  development  program  is  proposed  in  [41]:  an  extraction 
fluid  rate  of  140  kg/s  from  the  shallow  reservoir,  for  an  expected 
power  plant  with  maximum  size  3-3.5  MW  (back-pressure  or 
ORC),  serving  the  near  region  for  about  50  years. 


5.  Numerical  simulation  of  geothermal  reservoirs:  an 
“extended”  review  of  the  other  available  cases 

In  the  present  section  a  further  review  of  several  numerical 
simulations  is  reported.  The  cases  discussed  in  the  previous 
section  and  summarized  in  Table  1  deal  with  well  known  geother¬ 
mal  fields  and  concern  mainly  with  numerical  modeling  and  data 
matching.  The  numerical  simulations  considered  are  referred  to 
several  geothermal  fields  object  of  analysis  all  around  the  world 
[21,23]  and  [43-70],  The  additional  cases  listed  in  this  section 
often  deal  with  fields  which  are  presented  in  the  same  paper 
together  with  the  numerical  simulation.  Different  softwares  are 
used  and  various  grid  shapes  and  configurations  are  adopted. 

For  some  of  the  cases  reviewed  only  a  single  paper  is  available, 
being  the  experimental  data  on  the  field  not  available,  while  in  the 
most  recent  works  a  large  amount  of  experimental  and  explora¬ 
tion  data  are  available.  Different  geographic  areas  are  involved,  in 
order  to  perform  a  meaningful  worldwide  review.  Also,  different 
resource  types  are  considered  in  the  range  of  the  medium-high 
enthalpy.  The  cases  contained  in  this  section  have  been  reviewed 
in  the  Doctoral  Thesis  of  one  of  the  authors  (Vaccaro  [6]).  Table  2 
illustrates  the  main  characteristics  of  the  geothermal  field,  namely 
the  type  of  resource,  fluid  production,  plant  productivity  (only 
estimated  in  some  cases),  reservoir  extent  and  production  depth. 
Table  3  contains  the  main  indications  about  the  numerical  models, 
the  objectives  of  the  simulation,  the  geometry  and  extent  of  the 
model  and  grid  and  the  boundary  conditions.  Some  common 
aspects  can  be  found  in  the  boundary  conditions  and  calibration 
process  (due  to  the  uncertainty  about  permeability). 


6.  Discussion  about  the  numerical  models  reviewed  and 
perspective  for  utilization  in  the  energy  production  perspective 

As  it  is  evident  from  the  previous  review  analysis,  the  numer¬ 
ical  simulation  of  a  geothermal  reservoir  is  a  well  known  field  in 
the  literature.  In  Sections  4  and  5  a  summary  of  the  cases  is 
reported,  they  are  all  different  for  various  reasons.  First  of  all  for 
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Table  2 

“Extended  review”:  characteristics  of  the  geothermal  fields. 


Field  name  and  References 

Resource 

Production 

Plant :  Size 

Extent 

Heat 

Max. 

(kg/s) 

and  Type  (MW) 

(km2) 

flux 

prod. 

(mW  / 

Depth 

m2) 

(m) 

Balcova-Narlidere,  Turkey;  (Gok  et  al. 

Water  dominated,  140  °C 

20.6a 

0 

*2 

1100 

143]) 

(annual 

avg.) 

Cerro  Prieto,  Mexico  (Butler  et  al.  [44]) 

Water  dominated,  double  phase  reservoir,  350  °C 

620b 

3000 

Den  Haag,  Netherlands  (Mottaghy  et  al. 

Water  dominated,  single  phase  reservoir,  70  °C 

*42a 

5  (thermal,  estinated) 

547 

60-70 

2200 

[45]) 

Desert  Peak,  USA  (Zeng  et  al.  [21]) 

EGS  system,  *  210  °C 

*42 

2-5  (EGS-derived 
binary  plant,  under 
development) 

4.2 

(63) 

1200 

German  Basin  (Vogt  et  al.  [23]) 

Water  dominated  reservoir,  sandstone  reservoir,  87  °C 

%  42 

District  heating 

25 

50-78 

2000 

Hammam  Faraun,  Egypt  (Zaher  et  al. 

Water  dominated  reservoir,  70  °C 

Plant  for  desalinization 

77 

120 

[46]) 

Hammam  Musa,  Gulf  of  Suez,  Egypt 

Water  dominated  reservoir,  70  °C 

13  (estimated  potential, 

45 

80 

(Zaher  et  al.  [47]) 

volumetric  method) 

Hatchobaru,  Japan  (Yahara  and  Tokita 

Water  dominated,  double  phase  reservoir,  240-300  °C 

*700 

110 

16.5 

2300 

[48]) 

(*500 

reinjected) 

Heber,  USA  (Boardmann  et  al.  [49]) 

Water  dominated,  double  phase  reservoir,  200  °C 

120c  (52  effective) 

1830 

Hengill  Area,  Iceland  (Gunnarsonn  et  al. 

Water  dominated,  double  phase  reservoir,  300  °C 

330 

210 

2000 

[50]) 

Kamojang,  Indonesia  (Suryadarma  et  al. 

Vapor  dominated  reservoir,  230-245  °C 

300-424 

200 

14-21 

*1400 

[51]) 

Kizildere,  Turkey  (Yeltekin  et  al.  [52]  and 

Water  dominated,  double  phase  reservoir,  200  °C 

250  (10 

20.4 

1240 

Ozkaya,  [53]) 

Las  Tres  Virgenes,  Mexico  (Guerrero- 

Water  dominated,  double  phase  reservoir,  250-275  °C 

effective) 

10-14 

16.9- 

800- 

Martinez  and  Verma  [54]) 

20.3 

1000 

Los  Azufres,  Mexico  (Maldonado  et  al. 

Water  dominated,  double  phase  reservoir,  240-280  °C, 

*475 

188 

*20 

1500 

[55]) 

vapor  dominated  after  exploitation 

Mt.  Apo  -  Mindanao,  Philippines 

Water  dominated,  double  phase  reservoir,  300  °C 

104 

2500 

(Esberto  and  Sarmiento  [56], 
Emoricha  et  al.  [57]) 

Pauzhetsky-Kamchatka,  Russia 

Water  dominated,  double  phase  reservoir,  200  °C 

250  (38  to 

6.8 

*5 

63 

1000 

(Kiryukhin  et  al.  [58,59]) 

reinj.) 

Ogiri,  Japan  (Kunamoto  et  al.  [60], 

1st  reservoir  (400-200  m  asl):  water  dominated,  50- 

« 320 (243 

30 

432 

1800 

Itoi  et  al.  [61]) 

130  °C.  2nd  reservoir  (0  m  asl):  water  dominated,  double 
phase,  230  °C 

reinj.) 

Olkaria,  Kenya  (Ofwona,  [62]) 

Water  dominated,  double  phase  reservoir,  200-300  °C 

121 

*80 

1600 

Onikobe,  Japan  (Nakanishi  and  Iwai 

1st  reservoir,  water  dominated,  two-phase  (-200  to 

12.5 

*  1 

175 

1300 

163]) 

300  m  a.s.l.),  200  °C.  2nd  reservoir:  water  dominated 

(-1200  to  -700  m  asl),  250  C 

Poihipi  -  Wairakei,  New  Zealand 

Vapor  dominated 

55.5  (19.5 

55  (25  effective) 

(Zarrouk  et  al.  [64]) 

Sabalan,  Iran  (Bromley  et  al.  [65], 

Water  dominated,  double  phase  reservoir,  140-250  °C 

reinj.) 

50+50  (planned) 

*100 

*200 

*3200 

Strelbitskaya  and  Radmehr  [66], 
Noorollahi  et  al.  [67,68]) 

Sumikawa,  Japan  (Pritchett  et  al.  [69]) 

Water  dominated,  double  phase  reservoir,  200  °C 

0 

400 

2486 

Tanggu  District  (Guantao  field),  Tianjin, 

Water  dominated,  single  phase  reservoir,  60-75  °C 

11-25 

533 

*2000 

China  (Lei  and  Zhu  [70]) 


a  The  value  is  referred  to  a  district  heating  system  using  the  geothermal  resource. 
b  In  2000  an  upgrade  up  to  720  MW  was  planned,  a  further  upgrade  to  820  MW  was  planned  for  2012. 
c  Updated  at  1996  [49], 


the  kind  of  geothermal  field:  from  medium  enthalpy  water- 
dominant  field  to  dry-steam  dominant  field.  The  differences 
between  the  models  deal  with  simulation  domains  (size  and 
features),  scenarios  simulated  (unperturbed  or  exploitation)  and 
software  used. 

One  concept  has  to  be  emphasized  from  this  review:  the  strong 
dependence  of  the  results  of  the  numerical  analysis  on  the  quality 
of  the  inputs  and  the  difficulty  that  would  be  afforded  in  realizing 
the  models.  First  of  all,  the  data  and  the  geo-thermo-physical 
parameters  necessary  are  not  always  available  or  measurable. 
Moreover  defining  the  boundary  conditions  and  initial  conditions 
is  not  a  trivial  task. 

Initial  conditions  in  the  simulation  are  mostly:  current  thermal 
gradient  measured;  groundwater  recharge  flow  in  the  reservoir 


and  pressure  distribution.  The  boundary  conditions  deal  with  both 
the  hydro-geologic  and  the  thermal  problems,  regarding  essen¬ 
tially:  hydraulic  head;  temperature  distribution  at  the  top/bottom 
of  the  geometrical  domain;  heat  flow;  natural  recharge  of  steam  or 
water;  gases  diluted  in  the  geofluid  and  impermeable  and  adia¬ 
batic  conditions  at  lateral  faces.  All  of  these  conditions  can  vary 
with  time  during  the  simulated  interval. 

In  each  case  matching  between  the  forecast  and  the  measured 
data  is  then  fundamental.  The  measured  data  must  then  be  the 
basis  for  enhancing  and  improving  the  model  during  the  following 
steps  of  the  project.  For  most  of  the  listed  models  and  fields  it  is 
very  difficult  to  report  here  a  matching  between  modeled  reser¬ 
voir  and  real  measured  data.  In  many  cases  the  authors  of  the 
simulation  are  different  with  respect  to  the  utilization  plant  owner 


Table  3 

“Extended  review”:  numerical  model  simulations. 


Field  name  and 
References 

Simulation 

Software 

Geometry 

Mesh 

Conditions  and  parameters 

Calibration 

Balcova-Narlidere, 
Turkey  [43] 

Natural  state 

History  matching  (5  y) 
Evolution  (20  y,  3  cases) 

TOUGH2 

3D:  3  x  4  km2; 
thick.:  1.53  km  (30 
to  —1500  m  asl) 

Irregular  grid  13  layers  5194 
blocks3 

Deep  nat.  recharge:  40  kg/s;  shallow  recharge:  11  kg/s; 
total  heat  rate  input:  33  MW 

kb  natural  recharge 

Cerro  Prieto, 
Mexico  [44] 

Natural  state 

History  matching 
Evolution  (30  y,  3  cases) 

TETRAD 

3D:  11  x  10  km2 
thick.:  3.6  km 

7  layers  x  1296  cells;  9072 
blocks1;  range  size:  0.25- 
2  km 

Top:  7=40  °C;  Lateral  BC:  T  assigned  to  second  and  penultimate  layers, 
lateral  heat  flux,  fluid  losses  (347.2  kg/s);  Bottom:  T  assigned,  nat.  recharge 
347.2  kg/s  (350  °C);  </>  decrease  with  depth  (range  0.176-0.01),  fractures 
have  constant  porosity  kKy  assumed  as  a  logarithmic  function  of  $ 

k,  Xh  natural  recharge 

Den  Haag, 

Netherlands 
(Mottaghy 
et  al.  [45]) 

Natural  state 

Evolution  (50  y) 

SHEMAT 

ArgusOneIM 

(grid) 

3D:  22.5  x  24.3  km2 
thick.:  5  km 

2.43  x  106  nodes  9  geol.  units 

Top:  constant  T  (11  °C,  mean  annual)  Bottom:  constant  heat  flux  (63  mW/ 
m2);  X  depends  on  T 

Heat  flux,  X 

Desert  Peak,  USA 

Natural  state 

TOUGH2 

vertical  2D-3D; 

104  cells  ( x )  1  cell  (y)  104  cells 

Reference  case:  reinjection  1.5  kg/s  (60  °C);  <£= 0.2%;  k= 5  x  10-14  m2;  Mass 

Sensitivity  analysis  to  k,  X,  mass  flow 

(Zeng  et  al. 
[21]) 

Evolution  (20  y) 

height:  400  m; 
length:  400  m; 
thick.:  500  m 

(z)  (3  x  3,  3  x  5,  5  x  5) 

flow  rate  50-75  kg/s 

rate,  7rej,  po,  productivity  index 

German  Basin 

Natural  state 

SHEMAT-Suite 

3D:  5x5  km2  thick.: 

«  1.4  xlO6  cells  (18 

Top:  11  °C  surface  temperature;  Bottom:  regional  specific  heat  flow 

heat  flow,  temperatures  (3D  inversion); 

(Vogt  et  al. 

[23]) 

Evolution  (20  y) 

6  km 

sedimentary  layers) 

75  ±  10  mW/m2  (6000  m  depth)  Adiabatic  lateral  boundaries 

Stochastic  approach  to  hydraulic 
properties  assignment 

Hammam  Faraun, 

Natural  state 

HYDROTHERM 

3D:  77  km2  thick.: 

16  layers  x  (43  x  31 )  cells 

Top:  T,  p  constant  assigned;  Bottom:  heat  flux  fixed  (120  mW/m2);  Fracture 

Thermophysical  properties  determined 

Egypt (Zaher 
et  al.  [46]) 

(105  y+3  x  104y  for 
fracture  state) 

Evolution  (20  y) 

2.5  km 

(length  100-1000  m) 

zone  conditions 

by  previous  studies 

Hammam  Musa, 
Gulf  of  Suez, 
Egypt (Zaher 
et  al.  [47]) 

Natural  state 
(105  y+3  x  104y) 

HYDROTHERM 

3D:  45  km2;  thick.: 
3.6  km 

11  layers  x  (28  x  26)  cells 
(length  100-1000  m) 

Top:  T,  p  constant  assigned  (26.7  °C,  1.013  bars) 

Hatchobaru,  Japan 
(Yahara  and 
Tokita  [48]) 

Natural  state 

Evolution  (9  y,  3  cases) 

l 

3D:  16.5  km2;  thick.: 
2.5  km  (1100  to 
—  1400  m  asl) 

9  layers  (thick.:  100-400  m) 
7425  blocks 

Referred  to  a  trial-error  calibration  process  and 
further  papers  cited  in  [48] 

Trial-error  calibration  process 

Heber,  USA  [49] 

Natural  state 

TOUGH2 

3D:  14  x  13  km2; 
thick.:  3  km 

8  layers  x  201  cells  (1608 
blocks)*1 

Top:  7=25  °C,  p=0,l  MPa;  Lateral  BC:  impermeable,  adiabatic;  Bottom:  T 
and  p  assigned,  nat.  recharge 

kb 

Hengill  Area, 

Natural  state  (104y) 

TOUGH2 

3D:  100x100  km2; 

9  layers  x  996  cells  (8964 

Top:  7=15  °C,  p=0.1  MPae;  Penultimate  layer:  nat.  recharge  1  kg/s 

k  and  natural  recharge  (iTOUGH2) 

Iceland  [50] 

Kamojang, 
Indonesia 
(Suryadarma 
et  al.  [51]) 

History  matching  (20  y) 
Evolution  (30  y) 

Natural  state 

History  matching 
Evolution  (30  y) 

iTOUGH2 

thick.:  2.9  km; 

(400  m  to  -2500  m 
asl) 

3D:  49.5  km2; 
thick.:  3.6  km 

blocks)3 

15  layers  12480  blocks 

(1500  kj/kg)  Bottom:  7  and  p  assigned  (  as  265  °C),  fluid  flows 

Kizildere,  Turkey 

[52] 

Natural  state 

History  matching  (17  y) 
Evolution  (12  y,  9  cases) 

STARS  SAPHIR 

V.2.30 

3D:  840  x  600  km2; 
thick.:  1-1.2  km 
(estimated) 

5  layers  (8  x  12  cells) 

0,  k  (SAPHIR) 

Kizildere,  Turkey 

[53] 

Natural  state 

History  matching 
Evolution  (10  y,  4  cases) 

SUTRAf  v. 
1284-2D 

3D:  870  x  720  km2 

29  x  24  cells 

7  and  p  constant  at  boundaries  (not  specified) 

T.P 

Las  Tres  Virgenes, 
Mexico 
(Guerrero- 
Martinez  and 
Verma  [54]) 

Natural  state  (7  x  105  y) 
Emplacement  of 
magmatic  chambers 
Evolution  (30  y) 

TCHEMSYS 

3D:  20x30  km2; 
thick.:  22  km  (1900 
to  -20000  m  asl) 

I:  6307200;  blocks  II: 

4939200  blocks 

Boundary  temperature  at  the  value  given  by  the  gradient 
(0.03  °C/km),  surface  at  constant  25  °C 

Los  Azufres, 

Mexico  [55] 

Natural  state 

History  matching  (25  y) 
Evolution  (30  y,  2  cases) 

TETRAD 

3D:  12  x  12  km2; 
thick.:  3.2  km 
(3200-0  m  asl) 

9  layers  x  (18  x  26)  cells 
(4446  blocks )3,  c' s 

Top:  natural  emissions;  Bottom:  nat.  recharge  8.3  kg/s  (water  at  345  °C) 

Rock  parameters5  natural  recharge 

Mt.  Apo  - 
Mindanao, 
Philippines 

[56] 

Natural  state 

History  matching  (2  y) 
Evolution  (5  y,  2  cases) 

TETRAD 

3D:  6  x  10  km2; 
thick.:  2.75  km 
(1250  m  to  -1500  m 
asl) 

6  layers  x  (11  x  17)  cells  (1122 
blocks  )3,  h 

Natural  emissions  Bottom:  nat.  recharge  5  kg/s  (320  °C),  10  blocks 

k,  <t>  natural  recharge 
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Table  3  ( continued ) 


Field  name  and 
References 

Simulation 

Software 

Geometry 

Mesh 

Conditions  and  parameters 

Calibration 

Mt.  Apo  - 
Mindanao, 
Philippines 

[57] 

Natural  state 
(95  x  106  y) 

History  matching  (13  y) 
Evolution  (5  y,  2  cases) 

TOUGH2 

3D:  22  x  26  km2; 
thick.:  3.25  km 
(1250  m  to  -2000  m 
asl) 

19  layers  x  (31  x  47)  cells 
(27683  blocks,  only  16411 
active)3 

Top:  65  °C;  0.1  MPa;  Bottom:  nat.  recharge  145  kg/s  (320  °C);  Lateral  BC: 
impermeable,  adiabatic,  outflow  from  constant  pressure  cells 

k,  <p  natural  recharge 

Pauzhetsky- 

Natural  state 

TOUGH2 

3D:  polygonal 

131  cells  x  layer  Double  <p 

Top:  nat.  emiss.  100  °C;  Lateral  BC:  constant  T,  p;  Bottom:  heat  flux  63  mW/ 

k,  p  natural  recharge  rock  expansion 

Kamchatka, 
Russia  [58] 

History  matching  (36  y) 
Evolution  (20  y,  3  cases) 

A-Mesh 

thick.:  0.7  km  (avg.) 

model 

m2,  nat.  recharge11  204  kg/s  (830-875  kj/kg),  (SE)  120  kg/s  (900  kj/kg) 

coefficient 

Pauzhetsky- 

Natural  state 

TOUGH2 

3D:  polygonal1; 

3  layers  (424  blocks,  only  294 

Top:  atm.  T,  p;  nat.  emissions  (100  °C);  Lateral  BC:  impermeable,  constant  T, 

iTOUGH2:  T,  p,  natural  emissions. 

Kamchatka, 
Russia  [59] 

History  matching  (36  y) 

iTOUGH2 

A-Mesh 

thick.:  0.85  km  (100 
to  —750  m  asl) 

active)  Double  (f>  model 

p;  Bottom:  heat  flux  63  mW/m2,  nat.  recharge  224  kg/s 

natural  recharge  (m,  k)  History  matchng: 
production  wells n,  monitoring  wells  (T, 
p),  (p ,  k,  fractures 

Ogiri,  Japan  [60,61] 

Natural  state 

History  matching 

TOUGH2 

iTOUGH2 

3D:  5.5  x  3.9  km2; 
thick.:  2.85  km  (250 
to  -2600  m) 

7  layers  x  (23  x  11)  cells  (1771 
blocks)3, 1  cell  size: 

0.1-3  km  (thick.:  0.1 -1.6  km) 

Top:  T—  75  °C,  p=0.0981  MPa;  Lateral  BC:  impermeable,  adiabatic;  Bottom: 
heat  flux  43.2  mW/m2  (432  mW/m2,  S),  total  heat  in  19.5  MW  (260  mW/m2 
avg.);  nat.  recharge:6  30  kg/s  (240  °C),  total  inflow  (31.4  MW);  55  kg/s 
(1062.7  kj/kg,  inflow  58.4  MW),  production  area. 

k  (iTOUGH2) 

Olkaria,  Kenya  [62] 

Natural  state  (104y) 

TOUGH2 

3D:  polygonal 
(120  km2)m;  thick.: 
2.55  km  (2000  to 
—  550  m  asl) 

5  layers  x  158  cells  (790 
blocks) 

Top:  atm.  cond.e,  natural  emissions  (vapor)  366  kg/s;  Lateral  BC:  E-W 
impermeable,  N  p=45  bars,  S  p— 25  bars  (1075  m  asl);  Bottom:  nat. 
recharge  1253  kg/s  (6  blocks),  1600  kj/kg  (avg.);  fluid  loss  958  kg/s. 

kh  natural  recharge 

Onikobe,  Japan 

[63] 

Poihipi  -  Wairakei, 
New  Zealand 

[64,65] 

Natural  state 

History  matching  (21  y) 
Evolution  (10  y,  1  case) 
Prod.  Wells 

RANGER  STAR 

TOUGH2 

iTOUGH2 

AWTASn 

WELLSIM 

3D:  6.5  x  8  km2; 
thick.:  2.4  km  (400 
to  —2000  m  asl) 

14  layers  x  (10  x  11)  cells 
(1540  blocks)3,  gi  1 

Top:  const,  patm,  nat.  emiss.6;  Lateral  BC:  N-E  impermeable  and  adiabatic; 
S-0  constant  p;  Bottom:  heat  flux  175  mW/m2, 
nat.  recharge  10  kg/s  (330  °C)  production  area 

Porous/homogeneous  media  double  <p  (matrix/fracture)  fractional 
dimension  model 

kb  natural  recharge 

k,  cp  (AWTAS,  iTOUGH2) 

Sabalan,  Iran  [66- 
68] 

Natural  state 

Evolution  (3  cases) 

TOUGH2 

3D:  12  x  8  km2; 
thick.:  4.6  km 

16  layers  x  192  cells;  2688 
blocks  (2565  active)3,  b’  Cl  ° 

Top:  heat  loss  to  atmosphere;  nat.  manifestations  (modeled  as  4  artificial 
wells,  40-50  kg/s).  Second  to  last  layer  (near  to  the  base):  natural  recharge, 
8  kg/s  at  130  °C.  Lateral:  impermeable  and  adiabatic  Bottom:  nat.  recharge, 
90  kg/s  @  265  °C  (1159  kj/kg,  total  104.31  MW);  uniform  heat  flux: 

200  mW/m2  (Central-Southern  area),  null  in  the  North  area. 

K  natural  recharge 

Sumikawa,  Japan 

[69] 

Natural  state  (3  x  104  y) 
Evolution  (50  y,  1  case) 

STAR 

3D:  3  x  5  km2; 
thick.:  2.8  km  (1200 
to  - 1600  m  asl) 

16  layers  x  (9  x  10)  cells 
(1440  blocks)0,  p 

Top:  atm.  cond.  (T=10  °C),  nat.  Vapor  emissions  Lateral  BC:  E-W-S 
impermeable,  adiabatic,  nat.  recharge;  Bottom:  impermeable,  heat  flux 

400  mW/m2  (tot  6  MW)  Initial:  T  vertical  10-250  °C 

kh  recharge  geology 

Tanggu  District 
(Guantao 
field),  Tianjin, 
China  (Lei  and 
Zhu  [70]) 

Natural  state 

History  matching 
Evolution  (5  y,  2  cases) 

AUTOUGH2 

3D:  25x21.3  km2; 
thick.:  2.1  km 

2  layers  x  (780  blocks)  (thick.: 
150-1100  m);  Block  size: 
from  250  x  250  m2  to 

2000  x  2000  m2 

Water  table  assigned.  Top  layer  T  and  p  constant  (inactive  elements). 
Natural  recharge  and  inflows  are  assigned. 

P  thermophysical  parameters'3 

a  Mesh  refinement  in  the  production  wells  zone. 
b  Manual  iterative  calibration. 
c  Matrix  and  fracture  are  both  simulated  in  the  cells. 

d  Both  quadrangular  and  polygonal  shaped  blocks  (the  last  ones  are  used  to  simulate  the  fracture). 
e  Meteoric  water  recharge  is  also  considered. 
f  SUTRA  -  Saturated-Unsaturated  TRAnsport. 
g  First  layer  reproduces  orography. 
h  Mesh  refined  at  the  sea  level  layers. 

1  7  sides  polygonal  area:  2x0.5x3x2x3x  1.75  x  1.5  km7. 

1  M1NC  model  applied  to  the  main  fault  of  the  geothermal  system  (MINC  -  Multiple  INteracting  Continua). 
m  7  sides  polygonal  area:  11  x  11  x  5  x  9  x  10  x  2.5  x  6  km7. 
n  AWTAS  -Automated  Well-Test  Analysis  System. 

°  Mesh  refinement  in  the  atmosphere  contact  blocks. 
p  First  six  layers  reproduce  orography. 
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while  in  other  cases  the  history  matching  can  be  limited  to  a  quite 
reduced  time. 

Let  us  consider  two  different  cases,  the  Larderello  geothermal 
field  and  the  Sabalan  geothermal  field.  In  the  first  case  (see 
Table  1,  [30-33])  a  lot  of  historical  data  are  available,  also  from 
detailed  plant  production,  the  simulation  described  in  [30-32]  has 
been  realized  to  study  the  changes  in  the  large  field  after  years  of 
exploitation  and  drilling,  but  it  has  been  done  only  at  the  end  and 
not  in  a  forecast  mode.  For  this  very  famous  field  it  was  impossible 
to  think  about  a  strategy  as  the  one  proposed  here,  the  numerical 
simulation  being  a  modern  practical  solution  (the  first  plant  in 
Larderello  has  been  installed  in  1913).  In  the  case  of  Sabalan  (see 
Tables  2  and  3,  [65-68])  the  numerical  model  is  done  during  an 
exploration  step,  so  the  industrial  production  target  proposed  has 
not  been  verified  yet.  This  is  to  discuss  about  the  possibility,  for 
some  of  the  fields  and  models  illustrated,  of  matching  and 
comparing  simulated  and  measured  data. 

The  different  techniques  of  numerical  resolution  of  the  equa¬ 
tions  in  the  model  are  not  treated  as  a  problem  here.  But  it  is  clear 
that  it  is  a  common  problem  with  all  the  fields  in  which  numerical 
simulation  is  involved.  Although  the  codes  used  for  this  purpose 
are  evolving  very  quickly  and  the  results  can  be  very  detailed  and 
widely  complete,  these  simulations  present  a  remarkable  grade  of 
uncertainty.  In  the  perspective  of  a  more  diffused  industrial 
development  of  medium  to  low  temperature  geothermal  fields, 
the  numerical  simulation  can  be  a  very  useful  instrument  that 
must  be  connected  with  the  strategic  elaboration  and  environ¬ 
mental  and  economic  sustainability  of  the  design  of  a  geothermal 
plant.  Overall  a  relatively  good  agreement  was  obtained  in  the 
various  cases  between  the  measured  and  computed  temperature 
profiles,  but  the  matching  of  pressure  profiles  appeared  to  be  more 
difficult.  Even  if  the  definition  of  the  model  and  of  the  boundary 
conditions  requires  particular  attention  and  experience  in  order  to 
avoid  wrong  results,  numerical  simulation  could  be  a  good 
strategy  for  an  integrated  design  of  geothermal  plants  and  for 
the  prediction  of  the  geothermal  reservoir  response  and  environ¬ 
mental  impact  of  geothermal  plants  [71]. 


7.  Conclusions 

Numerical  simulation  is  a  fundamental  and  strongly  interacting 
instrument  for  plant  design.  Different  approaches  to  the  numerical 
simulation  of  geothermal  reservoirs  operation  are  considered  here, 
with  reference  to  some  case  studies  of  geothermal  fields  and 
ground.  The  perspectives  of  numerical  simulation  of  geothermal 
reservoirs  as  support  to  the  design  and  sizing  of  geothermal  plants 
are  outlined.  Models  simulation  is  a  powerful  decision-making 
tool:  it  can  provide  useful  indication  about  optimal  sizing  and 
sustainable  management  of  geothermal  utilization  systems. 

The  behavior  of  geothermal  reservoir  and  time  variations  of 
temperature,  hydraulic  head  and  pressure  have  to  be  estimated 
before  the  design  of  the  plant.  This  is  a  general  statement  for  each 
utilization  in  geothermal  energy  (steam  plants,  flash  plants,  dis¬ 
trict  heating  systems)  but  it  is  particularly  meaningful  in  case  of 
binary  plants  based  on  Organic  Rankine  Cycles  (ORC),  whose 
efficiencies  and  operations  strongly  vary  according  to  small  varia¬ 
tions  (concerning  mass  flow  rate,  temperature  and  pressure)  of  the 
characteristics  of  the  resource. 

Moreover  the  issues  of  scaling  and  of  correct  definition  of 
reinjection  strategy  must  be  considered.  The  importance  of  connect¬ 
ing  geological-geophysical  and  energy  engineering  background 
appears  fundamental  for  the  success  of  the  exploitation  of  a  geother¬ 
mal  field  in  order  to  sustain  the  geothermal  fluid  production  rate  over 
the  whole  lifecycle  of  the  plant  (in  general  higher  than  15  years). 


A  review  of  different  cases  of  numerical  simulation  from 
literature  is  considered  and  discussed.  The  state  of  the  art  methods 
and  commercial  softwares  available  for  the  simulation  of  geother¬ 
mal  fields  are  analyzed  through  the  review  of  several  geothermal 
fields  and  numerical  models  (ranging  from  medium  temperature 
geothermal  field  to  the  well  known  cases  of  Larderello  and  Mt. 
Amiata,  characterized  by  dry  steam  geothermal  resource).  Rela¬ 
tively  good  agreement  was  obtained  in  the  various  cases  between 
the  measured  and  computed  temperature  profiles;  simulation  of 
pressure  profiles  appears  to  be  more  difficult. 

Even  if  defining  the  model  and  the  boundary  conditions 
requires  particular  attention  and  experience  in  order  to  avoid 
wrong  results,  numerical  simulation  could  be  a  good  strategy  for 
an  integrated  design  of  geothermal  plants  and  for  the  prediction  of 
the  geothermal  reservoir  response  after  a  long  time  exploitation. 
Nevertheless,  the  following  issues  have  to  be  considered  in  order 
to  obtain  reliable  results. 

•  The  geothermal  phenomena  simulation  is  complex,  an  inter¬ 
disciplinary  approach  is  therefore  necessary. 

•  Software  and  codes  often  represent  solver  tools  for  the  physical 
equations  used. 

•  Case  study  experience  and  history  matching  is  a  fundamental 
background,  and  its  study  should  be  enhanced. 

•  Thermophysical  parameters,  boundary  conditions,  and  mesh¬ 
ing  method  strongly  affect  numerical  simulation.  Only  a  high 
accuracy  level  of  the  input  data  provides  reliable  results. 
According  to  the  quality  of  information  available,  simplified 
models  can  be  adopted. 

An  “integrated”  approach  to  the  complexity  of  the  geothermal 
phenomena  is  still  lacking.  Geothermal  energy  is  a  particular 
renewable  source:  its  use  is  sustainable  only  under  particular 
conditions,  which  must  be  known  particularly  by  investors  and 
market  players. 
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